Consider the orbit of a body moving under the influence of a spherically symmetric n-dimensional power potential with an arbitrary real exponent. The Hamiltonian for the system is

$H=\frac{p^2}{2m}+a{r}^k=E$

where the variables have their usual definitions as vector squares:

$r^2=(\mathbf{r}·\mathbf{r})=\sum _{i=1}^n{r}_i^2{p}^2=(\mathbf{p}·\mathbf{p})=\sum _{i=1}^n{p}_i^2$

Since the motion takes place in an invariant plane of two dimensions, the overall dimensionality of the space will not enter the development. The distinctive features of the system are not a result of our physical three-dimensional space but apply to any system with at least two spatial dimensions.

The straightforward approach to the integration of orbits is to replace momenta conjugate to angular variables with the square of angular momentum, thereby defining an energy equation for the radial variable:

$\begin{array}{r}\frac{1}{2}m{\stackrel{·}{r}}^2+\frac{L^2}{2m{r}^2}+a{r}^k=E\\ \sqrt{\frac{m}{2}}\int \frac{dr}{\sqrt{E-a{r}^k-\frac{L^2}{2m{r}^2}}}=t\end{array}$

Integration and inversion of this equation in principle produce the radial variable as a function of time. Angular momentum conservation

$m{r}^2\stackrel{·}{\phi }=L$

provides a differential equation linking a polar angular variable and the radial variable, integration of which produces the polar angular variable as a function of time.

The major obstacle to this approach is that the energy equation is rarely exactly soluble. It is notably so for the Kepler potential and the simple harmonic oscillator, which is why these two problems are worked to death in physics. The energy equation is also exactly soluble for a set of integral and rational exponents leading to elliptic functions or their combinations, but such solutions are not particularly transparent.

What one would like is a relatively simple statement of the general behavior of all orbits as a analytic function of the exponent of the power potential. One way to achieve this is with an approximate solution to a nonlinear differential equation. Although this is not an exact solution and also employs elliptic functions, it is relatively transparent and surprisingly accurate in the overall behavior.

Since the motion in the plane has only two spatial variables, label them the usual x and y for familiarity. Their equations of motion are

$\begin{array}{l}\stackrel{··}{x}=\frac{{\stackrel{·}{p}}_x}{m}=-\frac{1}{m}\frac{\partial H}{\partial x}=-\frac{ka}{m}{r}^{k-2}x\\ \stackrel{··}{y}=\frac{{\stackrel{·}{p}}_y}{m}=-\frac{1}{m}\frac{\partial H}{\partial y}=-\frac{ka}{m}{r}^{k-2}y\end{array}$

where again there is no loss of generality with respect to the n-dimensional space. A numerical integration of these equations for k > 0 with some typical initial conditions, where x is shown in red and y in green,

makes clear that the motion is always periodic for bound orbits, although the orbits are not in general closed. There appear to be two frequencies describing the motion: a fast frequency arising from the periodic motion of the radial variable, and a slower frequency describing the relative phase of the two spatial variables. This relative frequency becomes zero for k = 2 , the simple harmonic oscillator.

To capture this motion with the approximate solution to the nonlinear differential equation, first note from the energy equation that

$r\stackrel{·}{r}=\sqrt{\frac{2}{m}}\sqrt{E{r}^2-a{r}^{k+2}-\frac{L^2}{2m}}$

and evaluate the second temporal derivative of the square of the radial variable:

$\frac{d^2}{d{t}^2}{r}^2=2\frac{d}{dt}r\stackrel{·}{r}=\frac{2}{m}[2E-(k+2)a{r}^k]$

Scaling variables with the substitutions

$r^2=c_r^{2}gt=c_t\tau$

the second derivative becomes

$\frac{d^{2}g}{d{\tau }^2}=\frac{2c_t^2}{m{c}_r^2}[2E-(k+2)a{c}_r^k{g}^{k/2}]$

With the choices

$c_r=(\frac{2}{k+2}\frac{E}{a}{)}^{1/k}{c}_t=\sqrt{\frac{m}{4E}}{c}_r$

the statement of the second derivative is

$\frac{d^{2}g}{d{\tau }^2}+g^{k/2}=1$

which is the form given at the outset of the presentation of the nonlinear equation. The approximate solution from the presentation, with adjustments to variable names and an elliptic parameter rather than modulus, is

$g\left(\tau \right)\approx c_1+(c_0-c_1){cn}^2(\omega \tau ,\mathrm{m})$

with two of the three dependent parameters given by

$\begin{array}{l}{c}_1:\frac{c_1^{n+1}}{n+1}-c_1+l=0\\ \omega =\sqrt{2}K\left(\mathrm{m}\right)[\underset{c_1}{\overset{c_0}{\int }}\frac{dc}{\sqrt{c-\frac{c^{n+1}}{n+1}-l}}{]}^{-1},l=c_0-\frac{c_0^{n+1}}{n+1}\end{array}$

where $n=\frac{k}{2}$ and the elliptic parameter is denoted with an upright letter to distinguish it from the mass. The two constants c₀ and c₁ are extrema of the motion where $\stackrel{·}{r}=0$ . The first is related to an initial radial value by

$c_0=\frac{r_0^2}{c_r^2}=(\frac{k+2}{2}\frac{a}{E}{)}^{2/k}{r}_0^2$

This approximate solution always matches the exact radial variable at extrema. In the previous presentation, the form of the elliptic modulus was constructed to minimize the deviation of this approximation from a numerical solution between extrema. In this application a different choice will be necessary for the parameter.

To keep the solution to the nonlinear equation real, the first extremum is limited by

$1<c_0<(n+1{)}^{1/n}$

which in terms of an initial radial value and given energy is

$(\frac{2}{k+2}\frac{E}{a}{)}^{1/k}<r_0<(\frac{E}{a}{)}^{1/k}$

The lower limit represents a circular orbit at constant radius and the upper limit an orbit through the origin with zero angular momentum. These limiting orbits are described here in more detail.

For direct interpretation of orbits, the approximate solution needs to be written in terms of physical quantities. The square of the radial variable is

$r^2\approx r_1^2+(r_0^2-r_1^2){cn}^2(wt,\mathrm{m})$

where by definition $wt=\omega \tau$ . The quantity l appearing in the solution is

$l=(\frac{k+2}{2}\frac{a}{E}{)}^{2/k}\frac{L^2}{2mE}=\frac{1}{c_r^2}\frac{L^2}{2mE}$

so that the two dependent parameters are now

$\begin{array}{l}{r}_1:2mE{r}_1^2-2ma{r}_1^{k+2}-L^2=0\\ w=\frac{K\left(\mathrm{m}\right)}{m}[\underset{r_1}{\overset{r_0}{\int }}\frac{rdr}{\sqrt{2mE{r}^2-2ma{r}^{k+2}-L^2}}{]}^{-1}\end{array}$

Another input needed to capture the motion of the spatial variables x and y is the angular variable in the invariant plane. It can be evaluated from angular momentum conservation by a single quadrature:

$\phi =\frac{L}{m}\int \frac{dt}{r^2}=\frac{L}{mw}\int \frac{d\left(wt\right)}{r_1^2+(r_0^2-r_1^2){cn}^2(wt,\mathrm{m})}$

This integral can be expressed using a relation between the incomplete elliptic integral of the third kind and Jacobi elliptic functions:

$\mathrm{\Pi }[n;am(z,\mathrm{m}),\mathrm{m}]=\underset{0}{\overset{z}{\int }}\frac{du}{1-n{sn}^2(u,\mathrm{m})}=\underset{0}{\overset{z}{\int }}\frac{du}{1-n+n{cn}^2(u,\mathrm{m})}$

The approximation to the angular variable is thus

$\phi =\frac{L}{mw{r}_0^2}\mathrm{\Pi }[1-(\frac{r_1}{r_0}{)}^2;am(wt,\mathrm{m}),\mathrm{m}]$

The last item needed to capture the motion is a way to determine the elliptic parameter. Using the energy equation, the quadrature from the statement of angular momentum can be written in terms of the radial variable,

$\phi =\frac{L}{m}\int \frac{dt}{r^2}=\frac{L}{\sqrt{2m}}\int \frac{dr}{r^2\sqrt{E-a{r}^k-\frac{L^2}{2m{r}^2}}}$

so that one can define a change in the angular variable between the two extrema:

$\Delta \phi =\underset{r_1}{\overset{r_0}{\int }}\frac{Ldr}{r\sqrt{2mE{r}^2-2ma{r}^{k+2}-L^2}}$

This integral differs from the one determining the frequency of the radial motion in that the single power of the radial variable appears in the denominator rather than the numerator. This angular change corresponds to taking a complete elliptic integral in the approximate solution, i.e. setting its argument equal to $\frac{\pi }{2}$ . That means that if the elliptic modulus is chosen to satisfy the equality

$\frac{L}{mw{r}_0^2}\mathrm{\Pi }[1-(\frac{r_1}{r_0}{)}^2;\frac{\pi }{2},\mathrm{m}]=\Delta \phi$

then the approximation to the angular variable will always match the actual angular change at extrema, just as the approximation to the radial variable always matches the exact value at extrema. This is the best way to choose the elliptic modulus for this problem: although it allows error in the approximation between extrema, it always assures a match to the exact solution at the extrema.

Although one can write the elliptic modulus as an inverse complete elliptic integral of the third kind,

$\mathrm{m}={\mathrm{\Pi }}^{-1}[1-(\frac{r_1}{r_0}{)}^2;\frac{mw{r}_0^2\Delta \phi }{L}]$

in practice the value must be determined as a numerical root of the previous equality. While elliptic integrals of the first kind are inverted by the Jacobi amplitude function, elliptic integrals of the second and third kind currently have no convenient functional representation. Apparently even Wolfram is looking for such an inverse...

Note that for k < 2 the elliptic parameter satisfying this equality is negative and requires appropriate numerical handling, and Math is capable of that. This point will be relevant when considering k < 0 .

One can now present the complete approximate solution in an interactive graphic along with a numerical integration of the differential equation. Rather than restricting the radial variable for a given energy, the graphics will allow free variation of input parameters and calculate the energy from them. This is a more realistic exploration of the solution for a physical system.

The approximate solution with k > 0 for the spatial variable x along with the numerical integration is

and that for the spatial variable y is

In each case the numerical solution to the differential equation is plotted in the same color as at the outset, red for x and green for y, along with the approximation in blue. For a wide range of input values each plot appears as a single line, meaning the approximation is a nice fit to the numerical solution.

Since the angular variable in the invariant plane depends directly upon the frequency of the radial motion, it remains to understand the second frequency apparent in the graphic given at the outset. Unlike the radial variable, the angular variable has a continuous secular increase whose average rate can be determined from the time interval between extrema. The time between extrema, clearly the same for both radial and angular variables, is given by

$\Delta t=\underset{r_1}{\overset{r_0}{\int }}\frac{mrdr}{\sqrt{2mE{r}^2-2ma{r}^{k+2}-L^2}}$

If the angular change between extrema and the elliptic modulus have been determined, then solving the equality

$\frac{L}{mw{r}_0^2}\mathrm{\Pi }[1-(\frac{r_1}{r_0}{)}^2;am(w\Delta t,\mathrm{m}),\mathrm{m}]=\Delta \phi$

is equivalent to evaluating the integral but is duplicative and not as numerically stable. The value of the resulting amplitude is in general not equal to the $\frac{\pi }{2}$ used to determine the elliptic modulus, and this is the source of the ‘beats’ that appear between spatial coordinates when k ≠ 2 . Another way to express this is to note that the angular change between endpoints is also in general not equal to $\frac{\pi }{2}$ when k ≠ 2 : both are features of the same phenomenon.

If one converts the frequency of the radial motion to a quantity suitable for use with circular rather than elliptic functions,

$w_{\mathrm{radial}}=\frac{\pi }{2K\left(\mathrm{m}\right)}w=\frac{\pi }{2m}[\underset{r_1}{\overset{r_0}{\int }}\frac{rdr}{\sqrt{2mE{r}^2-2ma{r}^{k+2}-L^2}}{]}^{-1}=\frac{\pi }{2\Delta t}$

then the frequency of the beats is given by the difference

$w_{\mathrm{beats}}=w_{\mathrm{radial}}-\frac{\Delta \phi }{\Delta t}=\frac{1}{\Delta t}(\frac{\pi }{2}-\Delta \phi )$

To verify this expression, repeat the interactive graphic from the outset with the addition of sine and cosine envelopes at this frequency,

which match the numerical behavior quite nicely in terms of the beats frequency, although clearly the actual functions are more complicated than simple superpositions of circular functions. One can also evaluate the beats frequency as a function of exponent k for a restricted set of angular momenta

to verify that it does indeed become zero at k = 2 for all parameter values.

For negative values of the exponent, both the energy and coupling constant are negative for bound orbits. Since there is nothing in the approximate solution that depends essentially on particular values for parameters, one can use the same code as above with appropriate regions for exponent and coupling constant. Additional restrictions will be placed on all input parameters to ensure numerical stability of the output, but there is sufficient freedom of input to determine the viability of the method.

The approximate solution with k < 0 for the spatial variable x along with the numerical integration is

and that for the spatial variable y is

Variation of input parameters indicates that the fit is rather good for large angular momenta, but increasingly less exact for lower values. The approximation still matches the exact solution at extrema, as can be seen in this additional graphic,

but the approximation is clearly more qualitative than quantitative between extrema.

It is at first rather surprising that the approximate solution should continue to work at all for k < 0 until one remembers that the elliptic parameter begins to become negative for k < 2 . The region k < 0 simply continues this behavior, with larger and larger negative values necessary to satisfy the equality determining the approximation to the angular variable.

One might think to improve the fit of the approximation by allowing the radial and angular variables to have separate elliptic parameters. While this can help the approximation to the radial variable, it destroys the approximation to the angular variable. This latter approximation depends explicitly on the former and must share a single elliptic parameter determined by the complete elliptic integral in order to work at all.

Uploaded 2018.09.21 — Updated 2022.11.10 analyticphysics.com